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Abstract: We study the dynamical evolution of strongly coupled field theories, 
initially in thermal equilibrium, under the influence of an external driving force. We 
model the field theory holographically using classical gravitational dynamics in an 
asymptotically AdS spacetime. The system is driven by a source for a (composite) 
scalar operator. We focus on a scenario where the external source is periodic in 
time and chart out the response of several observables. We find an interesting phase 
structure in the response as a function of the amplitude of the source and driving 
frequency. Specifically the system transitions from a dissipation dominated phase, 
via a dynamical crossover to a highly resonant amplification phase. The diagnostics 
of these phases include the response of the operator in question, entropy production, 
energy fluctuations, and the temporal change of entanglement entropy for small 
subsystems. We comment on evidence for a potential phase transition in the energy 
fluctuations of the system. 
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1 Introduction and Conclusions 

The dynamics of quantum field theories driven far from equilibrium is a fascinating 
topic, owing to the complex interplay of quantum and statistical behaviours in the 
system. While a quantitative understanding of how field theories respond to non¬ 
linear external sources remains in general an open problem, in recent years one has 
gained some insight into such phenomena. 

On the one hand progress in this direction has been driven by experimental de¬ 
velopments which allow for a detailed study. For instance the ability to simulate 
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many-body dynamics in cold-atom systems has led to the opening of a new frontier 
in dynamical simulations, cf., [1] for a recent review. On the other hand, theoretical 
horizons have been broadened with the gauge/gravity duality providing an excel¬ 
lent arena to explore the dynamics of strongly interacting many-body systems using 
(classical) gravitational dynamics in a suitable limit (cf., [2] for a not so recent re¬ 
view) . Coupled with the development of excellent numerical algorithms for studying 
dynamical problems in AdS gravity [3, 4], the confluence of ideas and techniques 
provides an excellent opportunity to further our understanding of out-of-equilibrium 
dynamics. 

A much studied protocol in this context is the quantum quench dynamics, 
wherein one takes a system initially in equilibrium, typically in the ground state, 
and subjects it to external sources which change the subsequent dynamics by mod¬ 
ifying the Hamiltonian. The rate at which sources act on the system controls the 
features of the subsequent relaxation, assuming that the sources are non-vanishing 
for a finite amount of time. The analysis of such a quench protocol has benefited both 
from theoretical understanding using standard quantum field theory technology in 
low dimensions [5-7] and from a wide array of examples that have been studied holo¬ 
graphically in the recent past [8-21]. In most cases the interest is in the approach 
to equilibrium at late times and the rate at which various observables thermalize 
[16, 22-37]. Note that since we inject energy in the process of the quench, even 
an initially pure state will appear to be well approximated by a thermal ensemble 
asymptotically (assuming that the field theory dynamics are sufficiently ergodic). 

A slightly different but related scenario is one where we subject a system, again 
initially in an equilibrium configuration, to an external driving source which keeps 
doing work on it throughout the entire time period under study. More specifically, 
we will be interested in examining the behaviour when the initial state is chosen to 
be a thermal density matrix, so that one can simultaneously explore the response of 
a quantum dissipative system. For non-linear dynamical systems the response under 
such external driving can provide insight into the dynamics via the coherent build-up 
of the response. 

Classical analogs of what we have in mind are situations where we drive a 
(damped) pendulum steadily or subject a viscous fluid to external forcing. The 
latter is particularly apposite, for the problem we study can be thought of as a 
hot deconfined plasma of a planar gauge theory disturbed by an external source, as 
studied in the hydrodynamic context in [38]. Rather than letting the driving grow 
without bound, we will subject our plasma to a periodic driving by turning on the 
source for a relevant operator. One therefore has two relevant dimensionful param¬ 
eters characterizing the situation: (a) The amplitude of the external force whose 
scaling dimension is set by the conformal weight of the operator we exploit and (b) 
The frequency of the external driving. The third scale which is the temperature 
of the initial equilibrium state can be factored out, if we are interested in describ- 
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ing the dynamics for a conformally invariant system, which is most natural in the 
gauge/gravity context. This scenario was explored in [20], who carried out a per¬ 
turbative analysis for small amplitudes of the driving source. A related analysis of 
periodically driving a quantum system near a critical point was undertaken in [18] . 

Gravitationally the problem we study is the following: we have a Schwarzschild- 
AdS 4 black hole modeling our initial thermal density matrix of a three-dimensional 
CFT. At some instant of time on the boundary we turn on a periodic source for 
a relevant scalar operator, which we specifically choose to be of dimension 2 for 
simplicity^ The physics of the system is captured by examining the behaviour or 
various observables as we vary the amplitude A and the period P of the driving 
(measured e.g. in units of the initial temperature). We will in particular extend the 
perturbative analysis of [20] valid for A <C 1 to the non-perturbative regime A ^ 1 
for a wide range of driving frequencies. We find that the system naturally exhibits at 
least four different phases which are depicted in phase diagram Fig. 1; two of these 
(labeled Ilb and III) are non-perturbative in A. 

Before we describe the different phases, let us examine for a moment the physics 
of the gravitational system qualitatively. Initially we have a planar black hole in 
AdS 4 . When we turn on the scalar source, we are injecting energy into the bulk. 
This energy does work on the system and simultaneously heats it up. The latter is 
seen by the fact that some of the energy falls behind the horizon, which grows^ - 
this is the gravitational response to the disturbance of the plasma. However, in this 
process we also induce an expectation value for the operator whose source we tweak. 
When we disturb the system ‘slowly enough’, the operative parameter measuring this 
being the product of the amplitude and the period, the system has time to catch-up. 
This is the dissipation dominated regime indicated by I in Fig. 1. In this regime the 
injected energy falls behind the horizon with little fanfare. 

As we ramp up the disturbance, the plasma is driven more and more non-linear, 
with a dynamical cross-over visible as we move into phases Ila or lib of Fig. 1. Note 
that the entire non-linear dynamics in the system is induced by the non-linearities of 
gravity, for we model the system simply by a free (massive) scalar field. In this phase 
the response gets more in-phase with the source. It is amusing to contrast this with 
non-linear scalar dynamics; we find that in this phase we can model the scalar IPI 
effective action induced from the gravitational interactions to be well mimicked by a 
polynomial potential (see [17] for previous studies of self-interacting scalars in AdS). 

^ This choice turns out to have several advantages as the dual scalar being conformally cou¬ 
pled to gravity in the bulk allows a certain level of technical simplification in various holographic 
renormalizations we need to do to extract physical data. 

^ As we will be describing the dynamics of Einstein-scalar system with the scalar field satisfying 
the null energy condition, the areas of the event and apparent horizon (in the canonical foliation) 
have to grow monotonically - a consequence of the area theorem [39] (see [40] for an excellent 
overview). We will elaborate on this point in §2.3. 
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Fig. 1: The "phase diagram" of the driven holographic plasma characterized by the period 
(P) and amplitude (A) of the driving force, measured in units of the initial thermal scale Tq. 
There are four distinct regimes marked on the diagram which are explained in the main text. 
(Jin refers here to the imaginary (or in-phase) part of the conductivity defined in Eq. (3.3). As 
we move from southwest to northeast in the figure, the system is driven into a more non-linear 
regime: the crossing of the grey-dashed boundary is the turn on of the in-phase part of the 
conductivity din in regime II, and the crossing of the blue-dashed boundary signifies the entrance 
into the resonance phase of regime III i.e., \(/)'^^^\ oc. The character of the different regimes 
is further illustrated by displaying the phase portrait of the scalar operator (expectation value 
against source) used to drive the plasma. 


In this regime there is less dissipation; the entropy production by the growth of the 
horizon area is slowed down relative to region 1. The primary distinction between 
the two phases Ila and lib themselves is the lag in the response seen as the period 
is increased (hence the tilt in the phase portrait). 

For even larger disturbances, we enter region III, where the system response gets 
highly resonant and there is a steep growth in the response. As one might suspect 
this is the domain where the gravitational non-linearities are strongest and indeed 
one can check that such behaviour is not visible for a polynomially (self-) interacting 
scalar. In the course of our investigation we explore not just the phase portrait, 
but various other physical quantities of interest, such as the growth of entropy and 
dissipation in the system, the rate at which entanglement is produced, etc.. For 
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instance, region Ilb is characterized by enormous fluctuations in the energy of the 
system over a single period and continuous but non-differentiable behaviour in the 
entanglement entropy of a sub-system. 

Let us contrast our results with the analysis in the perturbative regime of small 
amplitudes undertaken in [ 20 ].^ As one can see from phase diagram Fig. 1 for small 
amplitudes, A <C 1 , one is largely in the dissipation dominated linear response regime. 
This is indeed consistent with the analysis of [20], who explore the dependence of 
observables on both the period of the driving as well as the dimension of the per¬ 
turbing operator A. As for us the latter remains frozen and we are unable to check 
the detailed scaling relations they And, but in the common domain of overlap we do 
indeed have agreement. In particular, for perturbing operators of dimension A = 2 in 
CFT 3 we expect to see that the energy dissipation as a function of the period scales 
as F'diss ~ P~^ (for A <C 1), independent of the initial temperature. Furthermore, 
we also expect that the work done in each cycle, measured by the entropy produced, 
to scale with the increased energy density. We And that in the slow driving regime 
this scaling closely tracks the prediction from local thermal equilibrium, but starts 
to grow more steeply as we transit into more interesting non-linear regimes. 

While the various response functions provide us with a useful diagnostic of the 
phase structure of the dynamical evolution, we also attempt to gain insight into the 
non-equilibrium dynamics using entanglement entropy for small sub-systems. This 
non-local probe exhibits distinct characteristic features in the various regimes: for 
weak driving, the growth of entanglement is gradual (and appears to track the growth 
of thermal entropy), while for strong driving there are steep oscillations and glitches 
in its evolution. We should caution the reader that we have only examined entan¬ 
glement entropy for relatively small sub-systems, owing to technical complications 
with numerical stability. Nevertheless these results suggest a rather rich structure in 
the temporal growth of entanglement with driving, which deserves further detailed 
exploration [41]. 

The outline of this paper is as follows. We begin in §2 by giving a quick overview 
of the basic set-up and the numerical solutions. Following this in §3, we set out the 
various observables we use to explore the behaviour of the system. In particular, we 
justify the rationale behind phase diagram Fig. 1 and how we should physically think 
of the different regimes. §4 is devoted to the study of entanglement entropy in this 
system where we focus on the region of an infinite strip and exploit the underlying 
homogeneity of the set-up. We conclude with a discussion in §5. Some technical 
results about holographic renormalization required for computing various observables 
is collected in the Appendices; Appendix A collects some useful information about 

^ We note that [18] study the influence of a periodic electric fields on the phase transition between 
a normal and superconducting phase using holography. It is clear in this case that a driving the 
system will make it exit the low temperature superconducting phase as the energy expended heats 
up the system past the critical point, as their analysis confirms. 
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holographic renormalization in our models while Appendix B provides details relevant 
for computing entanglement entropy. 

2 Driven CFTs and their Holographic Duals 

We first take the opportunity to set up the basic problem of a field theory driven 
out of equilibrium by turning on a source for a relevant operator. We then go 
on to describe how to model this in the holographic set-up and present the basic 
methodology and results from the numerical simulations. 

2.1 Driving CFTs by Relevant Operators 

We are interested in the dynamics of strongly coupled plasmas that are driven by an 
external source. The initial plasma is in equilibrium in some homogeneous thermal 
state at a temperature Tq for t < 0. At t = 0 we introduce external sources with some 
specified spatial-temporal profile that we control. We focus exclusively on situations 
where the external sources are spatially homogeneous, but otherwise arbitrary and 
tunable at will. 

To wit, the system under consideration can be modeled by an equilibrium density 
matrix, evolved under a time-dependent Hamiltonian, i.e., we take 

ScFT = Sj=o + j 7 J{x) 0{x) (2.1) 

where 0{x) is a single trace (gauge-invariant) relevant operator of conformal di¬ 
mension A < d. The source J7’(x) is chosen to have no spatial dependence and be 
temporally periodic and thus can be represented as 

J{x) = A cos{ut) 0(t). (2.2) 

Here 0(t) is the Heaviside step function for turning on the periodic perturbation 
of amplitude A and driving frequency (u = 27r/P at t = 0; later in actual (numer¬ 
ical) implementations we will choose a suitable ramp factor to smoothly turn the 
perturbation on. 

In the presence of the source, the Ward identities following from diffeomorphism 
and Weyl invariance get modified. A simple analysis shows that the boundary con¬ 
servation equation now has an explicit source term 

= (2.3) 

indicative of the work done by the driving source on the CFT. Likewise the one-point 
function of the trace of the energy-momentum tensor no longer vanishes but satisfies 

T; = (d - A) J{x)0{x) (2.4) 
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Since the boundary theory is conformal, it does not have any intrinsic time scale. The 
time scales in the problem come from only the driving force, namely its amplitude 
and period. The situation of interest is thus characterized by three scales: 

• Tq: the initial thermal scale for the homogeneous plasma. 

• A: the amplitude of the source whose scaling dimension is d — A. 

• oj: the driving frequency or the time-scale set by the period P = 27r/a;. 

2.2 Holographic Driving 

The gravity dual to this set-up is modeled by the dynamics of a scalar field (p with 
mass = — 2, dual to a relevant perturbation of the boundary theory. 

-S'buik = j [R + d{d-l)-^[ {dpf + (2.5) 

In our holographic implementation of this set-up we will work in d = 3 and consider 
a scalar operator with conformal dimension A = 2. While this is rather specific, we 
will explore the phase structure of the driven system as a function of the ratio of 
scales outlined above. The qualitative features we believe are independent of these 
actual choices.^ We have included a dimensionless gravity-scalar coupling ag which 
we can use to tune the amount of backreaction on the geometry; for the most part 
we will focus on ccj, = 0 or = 1, to model probe and interacting scalar fields 
respectively. 

We want to study gravitational dynamics driven by a scalar field whose non- 
normalizable mode is turned out as dictated by the source J{x), i.e., take p^it) = 
A cos{ut) and study the behaviour of the theory with varying amplitude A and 
frequency uj. The gravitational background is an asymptotically AdS 4 spacetime, 
which we write in ingoing Eddington-Finkelstein coordinates (sometimes called the 
Bondi-Sachs form) as: 

ds^ = -2 f{t, r) e^^d,r) de + 2 ^ ^^ 2 ^ ^ 2 . 6 ) 

The coordinate dependences of the metric functions ftXiP are explicitly indicated 
with homogeneity ensuring that dx and dy are Killing vector fields. 

Our initial state is a planar Schwarzschild-AdS 4 black hole with temperature 
To = S/tt, corresponding to horizon size r+ = 1. This bulk solution is given by 
/ = r^(l — ^), X = 0, p = r with metric 



r We have also set ^AdS = 1 for simplicity. 
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For our choices of = — 2 in d = 3, the amplitude A has mass dimension 1. Thus 
we have two interesting time scales associated with the external driving force: the 
period P and the inverse amplitude A~^. To capture universal physics, we look at 
relatively late times of the non-thermalized system compared to both of these scales. 
Note also that in those late times the initial value of the temperature, Tq, becomes 
irrelevant. 

There has been much interest recently in holographic quenches, in which the 
system is initially driven to an excited state, and then is allowed to return to equi¬ 
librium, a process which exhibits some degree of universality. In contrast, we are 
interested in the dynamics of the steady state system while it is being driven. Hence, 
in our solutions we do not turn off the driving force at late times, and seek univer¬ 
sal features associated with the driven steady state system. We will see that such 
dynamical features exist, and they strongly depend on the parameter 

i{P,A) = PA, (2.8) 

the unique dimensionless parameter formed from the two time scales associated with 
the driving force. Below we refer to the regime ^ <C 1 as the weak driving regime, 
and ^ 1 as the strong driving regime (which is further divided into two separate 
dynamical regimes). We also measure time in units of the period P, thus we vary 
and discuss the dependence of observables on the two dimensionless parameters: the 
strength of the drive and time. 

2.3 Bulk solutions 

We solve the equations of motion resulting from the scalar-gravity Lagrangian (2.5) 
by direct numerical integration. The boundary conditions on the scalar are pre¬ 
scribed by the source and the metric is required to be asymptotically AdS 4 . The 
AdS boundary is attained as r ^ oo and the asymptotic behaviour of the fields is 

r) — r + A(t) - ^ 

1 rv 

fit, r) = -(r + Xit)r - X'it) - -f Mtr + Oir-^) 

X(t,r) = C>(r"^) (2.9) 

More specifically, we use the characteristic formulation of the resulting partial differ¬ 
ential equations as explained in detail in [4] to numerically integrate for the solution. 
The advantage of the method is that it allows us to use constrained evolution: at 
each time step we solve a nested set of ODEs to determine the time derivatives of all 
dynamical quantities, and then we use one of the standard time evolution schemes 






(a) Sample </> solution (b) Sample / solution 


Fig. 2: A sample solution displaying the scalar field and the temporal component of 

the metric function f{t,u) for ^(P = 1,^ = 1) = 1. Time is measured in units of P and the 
radial component is compactified as u = 1/r. 


to march forward in time. While we follow the general logic of [4] , in our implemen¬ 
tation we found that some of elements described in [42] enabled for a more robust 
evolution. 

To solve the radial ODEs we discretize the equations using a Chebyshev basis in 
the radial direction, typically taking a grid of 60 points. For time evolution we use 
an explicit Runge-Kutta method of order 4, with an adaptive step size. We hlter at 
each time step by throwing out the top third of the Fourier modes for each dynamical 
variable to avoid artificial and unphysical growth in amplitudes of short wavelength 
modes associated with the UV cutoff. 

In the regime of strong driving, we found it necessary to turn on the perturbation 
gradually from zero. Therefore we include a ramp-up time of 2 P, after which the 
amplitude reaches its intended value. Thus, the first few periods of each solution 
show behaviour sensitive to details of the ramp-up protocol. We look at observables 
only after this ramp-up time of 2 P. 

In Fig. 2 we show one example of evolved bulk fields for a specific solution. As 
we perturb the system by a relevant operator, the scalar field grows towards the 
horizon. All fields are (at least approximately) modulated with the period of the 
source. 

At this point it is worthwhile mentioning one important consistency check on the 
numerical scheme, which relies on the existence of a smooth horizon in the spacetime. 
Given a metric and a Cauchy slice in the bulk spacetime, one can find the outermost 
trapped surface on this slice. If we have a set of Cauchy slices that foliate the 
spacetime, then the future outermost trapping horizon, which we simply refer to 
as the apparent horizon by a common abuse of terminology, is typically defined by 


9 








taking the union of the outermost trapped surface on all the slices. The apparent 
horizon thus defined is subject to an area law which was originally discussed by [39] - 
we refer the reader to [40] for a concise modern summary and proof of the statement. 
It is however important to note that the statement relies on the existence of a sensible 
foliation of the spacetime by Cauchy slices. Indeed, it is possible as discussed in [43] 
to find exotic symmetry breaking foliations (which are however incomplete) in which 
even the Schwarzschild black hole solutions fails to have a trapped surface. 

We mention this in passing, as [20] quotes the result of [43] to argue that apparent 
horizon areas need not be monotone generically. They however do not encounter such 
behaviour, for with the choice of ingoing coordinates in (2.6), there is a canonical 
choice of bulk Cauchy slices respecting the homogeneity of the disturbance. In this 
foliation the result quoted in [40] does apply and in fact simply follows from properties 
of null congruences using Raychaudhuri’s equation.^ Our results are indeed consistent 
with this expectation and we have checked that the area of the apparent horizon 
does grow monotonically in t (which labels the leaves of the foliation chosen), as we 
shall extensively see in the sequel. While initial results of [21] appeared to suggest 
otherwise, upon closer scrutiny, one finds that in numerical analyses so far the area 
of the apparent horizon does respect the second law as derived by [39].® 

3 Driving Diagnostics 

Having constructed the holographic duals we now turn to lessons that can be ex¬ 
tracted from the geometry for the dynamics of strongly coupled field theories. A- 
priori there are a number of observables which are useful probes of the out-of- 
equilibrium situation and we will focus on those that offer most clear insight into the 
dynamics. Our primary goal is to quantify the behaviour of the system as a function 
of {P, A} and construct a phase diagram demarcating the various regimes in this 
phase space. Let us quickly enumerate the observables we will use and proceed to 
explain why they give us some insight into the dynamics: 

• The phase portrait of response 0i(t) as a function of the source 0o(^)- Alter¬ 
natively, this relation can be codified in a conductivity cr(t), as defined below 
in (3.3). We find 4 underlying phases regions that the system can fall into. 

• The 01-00 phase portrait features for polynomial and non-polynomial potentials 
with the gravity-scalar coupling ag switched on and off. 

• The cycle-averaged thermodynamics quantified by the energy density eavg(t) 
and entropy density Savg(^)) and the scaling relation Savg ~ ^Zvg between them. 

® To be sure the statement of the area increase theorem does rely on the null energy condition, 
which we happily assume, for it is always satisfied by scalar fields with sensible kinetic terms. 

® We thank Alex Buchel for checking this and confirming the monotone growth of the apparent 
horizon area. 
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• The work done in each cycle, measured as the difference in average energy 
between two successive cycles, Ccyde = elvg^^ — elvg. We typically take n to 
correspond to the penultimate cycle of our simulation. 

• Fluctuations efluc(t) in the energy density around eavg(t) and the maximal re¬ 
sponse 

• Entanglement entropy and extremal surface evolution for fixed spatial strips A 
on the boundary. 

When the system is driven by an external source, the most basic quantity is the 
response, which is characterized by the scalar one-point function in the presence of 
the source. In linear response theory, this can be obtained from the retarded Green’s 
function of the operator 0{x) evaluated in equilibrium. We are not just interested 
in the linear response regime, which would correspond in our set-up to A -C Tq, but 
in the full non-linear response. To visualize the response of the strongly coupled 
plasma, especially in the non-linear regime, where its phase relative to the source 
is important, we will find it instructive to exhibit the phase portrait, the trajectory 
traced by the system in the 0o-0i plane. We also codify the relation between scalar 
source and response by a complex conductivity, defined below. 

In addition to the one-point function of the operator deforming the CFT, we are 
interested in the boundary energy-momentum tensor. This can be decomposed in to 
an energy density e{t) and a pressure. In the holographic set-up one has 

{Oit)) = Mt), e{t) = {T\{t)), pit) = {T\{t)) (3.1) 

The scale Ward identity (2.4) implies that pressure is not an independent observable 
since it can be obtained from knowledge of e{t) and 0i(t), so we will not discuss the 
pressure separately. Additionally, to probe the local thermodynamics we will monitor 
the local entropy density s{t), obtained by computing the area of the apparent horizon 
at time t7 

The dynamics of the bulk gravitational fields encode the heat production re¬ 
sulting from supplying external energy to the system. We monitor the explicit time 
dependence of the energy density e{t) and the entropy density s{t) along with their 
values averaged over each driving cycle period P, and find for the most part that the 
averaged values are increasing with time.® These provide a useful diagnostic of the 

^ Using the area of the apparent horizon (defined as the outermost trapped surface in the 
foliation respecting spatial homogeneity) results a causal boundary observable. One maps points 
on the apparent horizon to boundary points by Lie transport along radially ingoing null geodesics, 
which in the ansatz (2.6) are simply lines of constant {t,x,y}. On the other hand the teleological 
nature of the event horizon implies that its area would not provide a good measure for the boundary 
entropy density, cf., [3, 44] for a discussion of this point. 

® Note that the averaging makes eavg(t) and Savg(f) discrete in time. 
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departure from equilibrium, as one can monitor the scaling relation to infer the local 
thermodynamic equation of state. We define the thermodynamic scaling exponent 7 
when the system is in a steady state t > tg via 

Savg ^ ^avg (3-2) 

Note that in thermal equilibrium, conformal invariance predicts 70 = §• We will 
encounter this and other scaling regimes in our driven system when conformal in¬ 
variance is broken. 

Note that one natural set of non-local observables we could use are the multi¬ 
point correlation function for gauge invariant local operators, perhaps for O itself. 
However, realistically this computation involves solving the wave-equation for the 
linearized scalar fluctuations on top of the background we have constructed, to¬ 
gether with the imposition of suitable boundary conditions on the future horizon, to 
obtain sensible time-ordered correlation functions. These boundary conditions are 
somewhat tricky to implement (see however [45, 46]) - we will therefore postpone a 
discussion of correlators to the future.® 

Below we describe the behaviour of the observables mentioned above in three 
distinct dynamical regimes, and comment on the bulk interpretation of those regimes. 
Once we have gained sufficient intuition from this exercise, we will then examine the 
entanglement entropy for a specified boundary region. 

3.1 Dissipation Dominated Regime 

The simplest situation occurs in the regime of weak driving ^ -C 1, which is best 
described as the dissipation-dominated regime (phase I). This includes the regime of 
small amplitudes, studied perturbatively in [20]. In this weak driving regime, the 
behaviour of all observables is dominated by dissipation, which we now demonstrate 
by looking at some specific observables. 

As we drive the system by the scalar non-normalizable mode 0o it is instructive to 
divide the scalar response (pi to the part in-phase with the driving force, and the part 
completely out-of-phase with the perturbation. In analogy with an electromagnetic 
perturbation in linear response, we can complexify the time dependence of the scalar 
field^° and define a complex conductivity 

^ ^ crout(^) +icrin{t). (3.3) 

ICO cpoyt) 

^ We could following standard practice attempt to compute two-point correlation functions 
using the geodesic approximation [47]. However, as discussed in [48] and more recently in [49], 
this prescription doesn’t generically reproduce correct time-ordered correlation functions (we really 
want in-in correlation functions in our set-up). As a result we will also refrain from computing 
geodesics in the numerical background. 

That is, regard cos ujt and sincjt as the real and imaginary parts of 
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Fig. 3: The phase portrait of the dimensionless response ^ versus the dimensionless 
source 3 foi" ^ = 0.001 ^ 1 in the dissipation dominated regime {P = 0.001,^ = 1) 
which we label as phase I. We evolve the solution for 10 periods with each colour segment 
representing one period. The early times t <2P show the effect of the perturbation ramp-up, 
and thus are numerical artefacts that we omit from the plot. 


With this notation the out-of-phase and in-phase parts of the response correspond 
to the real and imaginary parts of the complex conductivity, crout(^) and crin(t), re¬ 
spectively. This is the usual convention for the more familiar conductivity, related 
to electromagnetic perturbations. As shown in Fig. 3 in the low driving regime the 
scalar response is precisely out of phase with the scalar source, (Tin = 0, meaning all 
the energy is dissipated and none of it used to excite the internal energy associated 
with the scalar field i.e., no work is being done on the system. This is the quench 
limit and it matches with what we expect from the behaviour of the perturbation in 
linear response. The complex conductivity a = (Tout is purely real and has constant 
amplitude as a function of time at high frequencies.^^ This is manifested in the fi¬ 
nal steady state being reached almost immediately and consisting of closed untilted 
trajectories in phase space. As we shall see below, tilting of the trajectories in phase 
space is indicative of non-trivial response and work done onto the system. Fig. 4 
shows what fraction of the complex conductivity (Tout is present on each point on the 
(P, A) phase diagram, and for what we are concerned with currently, the system has 
the response being completely out-of-phase with the source when the period is low. 

Both the energy and entropy density, averaged over each cycle, grow linearly 
with time in the dissipation-dominated regime . As the black hole grows, its entropy 
growth tracks its energy growth at a slightly higher rate than the equilibrium relation 
•Savg Cavg, i-G., 7 > 2/3. TMs entropy-energy scaling is shown in Fig. 5 along with 

This is similar to the behaviour of the conductivity for electromagnetic perturbations in asymp¬ 
totically AdS space. 
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Fig. 4: The fraction of the complex conductivity a\n over the entire (P^A) phase diagram 
where \a\^ = cr? + 



Fig. 5: The fitted average entropy Savg versus the average energy eavg (left) and their individual 
values as a function of time (right) for ^{P = 0.01,^ = 1) = 0.01. Fitting for Savg ~ eLg. we 
find a fitted value of 7 = 0.6682 zb 0.0023 > | with 95% confidence. 


their own evolution with time. Note that the expansion of the black hole horizon is 
not necessarily adiabatic (as measured e.g., by the rate of entropy increase |;|). 

In the low amplitude regime, one can also estimate in perturbation theory the 
amount of energy dissipated per cycle Ccyde which we dehne as the difference of the 
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1.6 



Fig. 6: The dimensionless scaling parameter a{Lo) from fitting Ccyde ~ foi" a small amplitude 
A = 1 m the linear response regime. It is expected for our choice of the scalar and dimension 
(A = 2 and d — 3) that cr —> 1 in both the small (ccycie ~ and large frequency (ccyde ~ 
limits. 


average energy e^vg between two successive cycles; for simplicity we take the result 
for the last two cycles of our evolution in quoting the results below. One expects the 
relation to take a scaling form Ccyde The scaling exponent a should be a non¬ 

trivial function of frequency itself; for low frequencies it is independent of the driving 
operator, but the high frequency limit cares about the spectral properties about the 
operator in question. Specihcally, one finds that [20] : Ccyde ~ for small frequencies 
and Ccyde ~ for high frequencies. Since we are not scanning over different 

choices of the driving operator, we have a single shot at determining this result. As 
depicted in Fig. 6 we indeed find that the energy dissipated is linear both at low and 
high frequencies: a{u)) —)■ 1 both for a; ^ 1 and for a; <C 1 (a coincidence owing to 
our choice A = 2 and d = 3). Interestingly there is some non-trivial intermediate 
frequency behaviour which appears to amplify the energy dissipated in a single cycle. 

The bulk picture of the process is also very simple: as we send energy pulses, 
which are either weak or infrequent, they interact very rarely before falling into the 
black hole horizon. All injected energy from the boundary goes towards steadily 
increasing the black hole mass and the scalar field remains unexcited. The more 
diverse behaviour observed below can be attributed to gravitational interactions of 
those energy pulses before they fall into the black hole. 

3.2 Dynamical Crossover Tilted Regime 

We now discuss the qualitative changes in the system as we begin to move from the 
weak driving ^ <C 1 to the strong driving regime ^ ^ 1 (from regime I to regime II 
through the grey-dashed line in phase diagram Fig. 1). Fig. 7 depicts a typical phase 
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portrait of the system as we cross into the new dynamical regime. We see that this 
regime is characterized by an onset of excitations of the scalar field and breaking 
of discrete time translation symmetry. The left panel of Fig. 7 shows the transition 



(a) ^(P = 0.1,T = 20) = 2 


(b) ^(P = 10,T = 1) = 10 
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Fig. 7: The dimensionless phase portrait of the response versus the source 0o ^{P = 
0 .1, A = 20) = 2 (left) and ^{P = 10, A = 1) = 10 (right). The conventions are as in Fig. 3. 
The left panel shows the behaviour in phase Mb while the right panel pertains to phase lla. 


from (^<Cl^<^^lat high amplitudes: the trajectories are no longer closed, rather 
they precess as a function of time and are slightly tilted. The breaking of discrete 
time-translation invariance is an interesting effect of the gravitational interactions of 
the scalar field. 

In the right panel of Fig. 7 we see the effect of moving into the new dynamical 
regime at low amplitudes: there is a clear tilt in the phase portrait from the one in 
Fig. 3 with ^ ^ 1 which indicates that the response is no longer completely out of 
phase with the source. The tilting of the trajectories at lower frequencies corresponds 
to the emergence of a finite in-phase contribution din > 0 in the conductivity; this sets 
the system somewhere between one with a purely out-of-phase conductivity (closed 
circular trajectories) and one with a purely in-phase conductivity (straight diagonal 
line trajectories). In other words not all of the injected energy is dissipated as was 
the case in regime I, but rather, work is actually being done on the system. 

As a result of having less dissipation in this regime, the energy and entropy of 
the black hole grow more slowly with time. Moreover, we find the scaling behaviour 
between the average energy and entropy, with a thermodynamic scaling exponent 
7 > |, for all values of A), as shown in Fig. 8. In other words, while the work 
done in the system slows down the energy increase of the black hole, the entropy 
production is affected less. 
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Fig. 8: The increase in the scaling exponent 7 in Savg ~ ^avg from the equilibrium value of 
7 o = I over the entire {P,A) phase diagram. We find that 7 > 70 holds for all scanned values 
on the phase diagram. 


To understand this regime further, it is instructive to reproduce this type of 
phase portrait for a system without gravity. To that effect, we can study the special 
case of scalar field evolution in a fixed black hole background, with no backreaction 
on the geometry (i.e., = 0). To include non-linearity into the problem, we add 

self-coupling to the scalar field, to mimic the effect of the non-linearities due to 
gravitational interactions (see also [17]). Fig. 9 depicts the phase portrait of a self- 
coupled scalar field with two types of polynomial potentials, which we took to be our 
original form (free massive scalar) and also one with quartic self-interactions: 

W(</>) = -20' - (3.4) 

We can see that without non-linearity as in Fig. 9a, the phase portrait is tilted, 
but sharp features of the phase portrait are lost compared to the case with the same 
driving but also gravitational backreaction, depicted in Fig. 7b. Adding a polynomial 
non-linearity, as done in Fig. 9b, gives a phase portrait that starts to form slightly 
sharper features along with some amplification of the response. Thus, the simple 
system of self-interacting scalar field allows us sufficiently separate the two effects 
in regime II: we see that the tilt in the phase diagram is associated with decreased 
frequency, whereas the breaking of time-translation invariance is associated with 
increased amplitude. We note also that for this simple system, the third dynamical 
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Fig. 9: The phase portrait of the dimensionless response ^ versus the dimensionless 
source 0 o = 3 for ^{P = 10 , ^ = 1 ) = 10 with ag = 0 and different polynomial potentials 
V{4>). The conventions are as described in Fig. 3. 


regime of unbounded amplification discussed in the next subsection seems to be 
absent. 

Thus, the bulk interpretation of this dynamical regime becomes clear: the pulses 
of energy injected at the boundary interact gravitationally before falling into the 
black hole. This results in additional physics to that of simple dissipation, modeled 
here by infalling the black hole. The gravitational interaction is due to perturba¬ 
tive exchange of gravitons, and can be mimicked by a polynomial self-interaction of 
the scalar field. In the next subsection we will see the effect of the gravitational 
interactions becoming strong when both A and P are large. 

3.3 Unbounded Amplification Regime 

As we increase the driving strength further in both A and P directions (from regime 
II to regime III through the blue-dashed line in phase diagram Fig. 1), we enter a 
dynamical regime no longer reproducible by polynomial self-interactions of the scalar 
field. We see the phase portrait of the scalar field in Fig. 10 for two instances of pa¬ 
rameters in this regime. Moreover, we find this dynamical regime to be characterized 
by unbounded response and restoration of time translation symmetry. 

As we increase the strength of the driving force the phase portrait becomes 
sharper and tilted, corresponding to an increased response and, again, less lag with 
the source as seen in Fig. 4. The ‘slowness’ of the energy injection from the boundary 
allows the scalar field to heat up as if the entire process were adiabatic, consequently 
allowing the scalar response to respond relatively quicker to the source. Note that 
although Fig. 4 shows |(Tin/cr| fa 1 in this regime, the absolute value |(t| is actually 
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(a) ^(P = 1,^ = 20) = 20 


(b) ^(P = 10, ^ = 20) = 200 


Fig. 10: The phase portrait of the response (pi versus the source (po for ^ = 20 (left) and 
^ = 200 (right) in the non-perturbative dynamical regime (regime III). The conventions are as 
in Fig. 3. 


very large in this unbounded amplification regime so that a small IcTout/c^l is still 
strong enough to keep the black hole perpetually growing in size. 

The maximal response over our ten cycles of driving is plotted in Fig. 11 

throughout the phase diagram. It is seen to increase rapidly with ^ past the dissipation- 
dominated regime. This seems to indicate the presence of a non-linear resonance, 
which allows the scalar response to grow without bound. An interesting feature of 
Fig. 11 is that the maximal response does not grow in the high frequency regime 
regardless of how large ^ is by increasing A. It seems unlikely that unbounded 
behaviour is attainable even for amplitudes drastically higher than the bounds of 
numerical explorations reported in Fig. 11. Physically, this means that a rapid puls¬ 
ing of small packets of energies can barely amplify the response of the system; the 
frequency of driving has to be below a certain bound for resonance to be possible 
- or in other words, a certain slowness in the sourcing is required. We conjecture 
that one should would see unbounded amplification only in the combined large P, 
large A regime which is slightly different from the traditional definition of resonance 
that depends only on frequency. An interesting curiousity is a slight dip in the re¬ 
sponse for moderate values of ^ preceding the rapid growth. This trough appears 
to demarcate the domains of bounded (regime II) and unbounded responses (regime 
III) empirically. It would be interesting to come up with a explanation for this 
phenomenon. 

Finally, it is amusing to model the non-linear effects of gravity in terms of an 
effective scalar potential to see what is necessary to attain regime III. We find that 
while a scalar field with polynomial self-interaction does not seem to posses this 
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regime, one can reproduce similar features by non-polynomial potentials. For exam¬ 
ple, we can discuss a self-interacting scalar field probe, with 

Kon-poiy(0) = -2 sinh^ 0 -F ^ sinh^ 0 . (3.5) 

This choice of scalar self-interaction is chosen to agree with our previous example 
(3.4) in the small field regime, but of course behaves differently for large field values. 
In Fig. 12 we see that indeed similar features of the phase diagram are reproduced: 
narrow closed trajectories and resonant response. We conclude therefore that the 
features of this dynamical regime are due to strong, non-perturbative gravitational 
effects occurring outside the black hole horizon. The fact that the non-linearities 
induced by gravity can be extremely strong, should perhaps be borne in mind while 
attempting to come up with simplified models of gravitational dynamics in AdS 
spacetime. 

3.4 Energy Fluctuations 

Another observable we monitor is the behaviour of energy fluctuations. More pre¬ 
cisely, we consider the deviations from the average energy in a each cycle, efluc(t) = 
|e(t) — eavg(t)|- These cycle fluctuations are a crude proxy for genuine fluctuation 
information that can be extracted, for instance, by considering symmetrized two- 
point functions of the boundary energy momentum tensor. Such ensemble-averaged 
fluctuations are known to exhibit phase transitions in periodically driven systems 
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(a) ^{P = 10,A = 1.5) = 15 


(b) ^(P = 10,^ = 2) = 20 


Fig. 12: The phase portrait of the response versus the source (po for ^ = 15 (left) and 
^ = 20 (right) for the non-polynomial potential Eq. (3.5), in the conventions of Fig. 3. 


[50]. Some indication those transitions are possible in holographic systems is given 
in [20]. 

The results for our simulations in various regimes are plotted in Fig. 13. We 
observe a qualitative change in these cycle fluctuations between different regimes. 
While in the dissipation-dominated phase we do not see a lot of deviation from the 
mean, there is a steep growth in fluctuations as we enter the non-linear phases. The 
fluctuations are maximal in the unbounded amplification regime (regime III). We 
note that in contrast to the maximal scalar response, which also grows dramatically 
in that phase, the fluctuations do track the driving frequency, with there being more 
deviations in the large period limit. 

It would be useful to confirm this behaviour directly with the computation of 
correlation functions, a task we leave for future investigation. 


4 Entanglement entropy 

Thus far we have discussed various local observables (response functions and ther¬ 
modynamic data) which have served to help us chart the phase diagram of the driven 
system in Fig. 1. We now turn to other non-local held theory observables that are 
sensitive to the non-equilibrium dynamics. Since we are not going to examine the 
behaviour of higher point correlation functions, we will dive right into the dynamical 
behaviour of entanglement entropy. 

In the boundary we have a density matrix p{t) which is time-evolving with respect 
to the perturbed Hamiltonian. At any given instant of (boundary) time, we pick a 
spatial region A and construct the matrix elements of the reduced density matrix 
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Fig. 13: Energy density fluctuations (last cycle) Cfiuc in units of /P over the entire (P, 
phase diagram. 


PA{t) = Tr^c {pit)) by tracing out the degrees of freedom in the complement (on 
the chosen Cauchy slice). The entanglement entropy is given by the von Neumann 
entropy of i.e., <S'^(t) = —Tr^ (p _4 logp^) which we can monitor as a function of 
time. 

Holographically computing the entanglement entropy for boundary regions in 
time dependent situations involves hnding bulk codimension-2 extremal surfaces 
anchored on the said boundary region A [23] . We study the evolution of entanglement 
entropy focusing in particular on translationally invariant strip regions: 

A = {t = —a < X < o, p e M} . (4.1) 


The bulk codimension-2 surface ends at x = ±a at some chosen instance of boundary 
time and is obtained by solving effectively a set of geodesic-like equations with 
our interpolated metric functions S, /, and y (see Appendix B.l for details). The 
covariant holographic entanglement entropy prescription [23] generalizing [51, 52] 
states that 


^ _ Area(£:^) 


(4.2) 


Should there be multiple extremal surfaces, we choose the one with minimal area 
(homologous to A). The proper area of these surfaces diverges owing to the locality 
of the underlying QFT. In our case we encounter potential divergences not only from 
the surface reaching out to the asymptotic boundary, but also from the presence 
of the sources driving the system. The physical result we are after is the finite 
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universal contribution which will measure the entanglement created/destroyed 
as we drive the system away from thermal equilibrium. Fortuitously, for our choice 
of scalar operator, there are no contributions due to the source, and hence we can 
simply regulate by background subtraction.^^ As a result we will consider as our 
entanglement diagnostic, the following finite quantity 

4 

[SJt) - SJt = 0)] (4.3) 

Liy 

where Ly is the IR regulator in the non-compact translationally invariant direction. 
Since we drive the system away from thermal equilibrium, = 0) is the corre¬ 
sponding value of the entanglement entropy computed in the Schwarzschild-AdS 4 
geometry. In what follows we will simply quote the results of our numerical simula¬ 
tions both for the behaviour of the extremal surfaces themselves and AS^ {t). 

4.1 Extremal surfaces in the driven geometries 

The extent to which the extremal surfaces penetrate into the bulk can for the most 
part be determined from the location of the cap-off point which we parameterize as 
{t*,u* = l/r*,x = 0).^^ For very small regions we are reasonably close to the AdS 
boundary whence, the curves are approximately semi-circles ~ a^. As we 

increase to larger strip widths the extremal surfaces start to probe the interesting re¬ 
gions of the driven geometry and thus allows us to see qualitative differences between 
the four phases. 

Generically we see that the following statements hold irrespective of the phases 
we consider: 

1. The radial depth and the temporal extent spanned by the surface evolves non- 
trivially as a function of One consequence of working with ingoing coordi¬ 
nates (2.6) is that the surfaces naturally dip back in time (see [33, 53]). 

2 . The oscillatory driving of the system imprints itself in the profile of the extremal 
surfaces, with the scale of these oscillations set by the the driving parameters 
A and P. The periodic movement of the surface can be seen in pulsations of 
the turnaround point of the surface: u* and t* have oscillations of the same 
period superposed over some enveloping function. 

3. On average, the extremal surfaces reach further into the bulk with time; 

is monotonically increasing for the range of parameters explored. To under¬ 
stand this note, we gauge fixed the bulk coordinate chart (2.6) such that the 

Details of the divergent structure and the counter-terms necessary to compute the area func¬ 
tional in our set up can be found in Appendix B.2. 

The coordinate u = 1/r is chosen such that the horizon remains at u = 1 during the entire 
course of the evolution (the boundary is at u = 0). 
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horizon is at = 1. In these coordinates the proper size of the region A 
increases (due to p(t,r)) which means that the surfaces want to get closer to 
the horizon to extremize the area functional. The rate at which this happens 
depends on both the amplitude and the frequency of the driving. We also note 
that surfaces dip less temporally, i.e., t* — is increasing. 

4. We also note that the location of the extremal surface appears to be consistent 
with causality of entanglement entropy [49]. While we have not explicitly 
checked that the surface lies in the casual shadow of the boundary region A, 
one simple consistency check visible from our results for t* is that t* — a. 
We remind the reader that in (2.6) lines of constant t and x are radially ingoing 
null geodesics. Causality at the very least requires that the cap-off point of 
the extremal surface lies below the ingoing null geodesic from the domain of 
dependence. Since for the strip region the boundary domain of dependence is 
a diamond anchored at (t^ ± a, 0) and (0, ±a), we note that the ingoing light 
ray from the bottom tip of this diamond cannot signal to the cap-off point. 

In the following discussion we will illustrate the behaviour of the extremal sur¬ 
faces more explicitly in each of our phases. We have been reasonably conservative 
in our analysis and have chosen to work only with surfaces that do not get too close 
to the horizon (in fact u* < 0.2). This is to avoid both numerical issues as well as 
to avoid complications from the existence of multiple extremal surfaces. We follow 
a single branch of solutions as described at the end of Appendix B.l. The primary 
results of the extremal surfaces are shown in the plots Figs. 14, 15, 16, and 17, where 
we show the evolution of the extremal surface as well as and 

Linear regime (small A): Although all phases display extremal surfaces that sink 
into the bulk with each driving cycle, the growth of u* in the linear regime of small 
amplitudes is most steady. We focus here on phases I (high frequency; dissipation- 
dominated) illustrated in Fig. 14 and Ila (low frequency; tilted) illustrated in Fig. 15, 
which fall under this characterization. As the frequency is lowered and we pass from 
the dissipation-dominated phase to the tilted phase, there is drastic reduction in the 
growth of u* per cycle. 

The evolution of t* in the two phases is also interesting; t* — is gradually 
increasing on average with time (recall that in the stationary geometry t* — would 
be constant). It turns out to be useful to look at a dimensionless parameter i* = 
(t* — which measures the cap-off time relative to the boundary. In this context, 

there is more time-lag in phase I i.e., t* <C < 0, which hints at the cause for why 
the surfaces do not penetrate as far deep in the bulk in phase Ila as opposed to phase 
I.^^ In addition we see strong oscillatory patterns in phase II in spite of having only 
a steady increase in u*; such a feature is absent in phase I. 

^^Note that in absolute terms however, t* in both regimes is comparable in magnitude. 
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Fig. 14: Evolution of the extrennal surfaces for a strip of width a = 0.05 with driving parameters 
^{P = 0.1,^ = 1) = 0.1 (phase I; dissipation-dominated). We pick a UV cutoff = 10“^ 
and have defined i* = {t* — tX)/P to measure the cap-off t* point relative to the boundary. 


Non-linear regime (large ^4): We now turn to the phases III (low freqnency; 
nnbounded amplification) illustrated in Fig. 16 and Ilb (high frequency; wobbly) 
illustrated in Fig. 17 in the non-linear regime of high amplitude. Some of the features 
seen in the linear regime continue to pertain: we see more pronounced oscillations in 
t* and a decreased tendency for the surfaces to lag behind in time at lower frequencies. 

In the unbounded amplification regime (phase III), we see significant bursts of 
growth of the extremal surfaces. The oscillatory driving is felt rather acutely by the 
surfaces and the evolution is considerably violent. On average however, u* appears 
to advance more serenely despite having large amplitnde oscillations per cycle. 

In the dynamical crossover wobbly regime (phase lib), there is a considerable 
amount of instability. We chose here to work with smaller strip widths a = 0.01 
(instead of o = 0.05) to avoid complications of phase transitions between multiple 
competing extremal snrfaces. The early part of the evolntion is in line with what 
happens in the dissipation-dominated regime (phase I), but shortly after, there are 
discontinuities in the i* parameter with no noticeable effect in u*. Around ^ 

4.0 — 4.2 and t^/P ~ 4.6, we see an exchange of dominance in the extremal surface. 
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Fig. 15: Evolution of the extrennal surfaces for a strip of width a = 0.05 with driving parameters 
^{P = 10,^ = 1) = 10 (phase II; tilted). Conventions described in Fig. 14 apply. 


which starts out at a higher value of t*. 

All in all, the extremal surfaces in the non-linear regime definitely has elements of 
intrigue owing to the large pulses of energy that affect the bulk geometry significantly. 
Although we do not delve into extremal surfaces that are positioned deeper into the 
bulk, we notice in the course of our analysis that the surfaces tend towards the 
horizon as expected. More curiously, we also find that for larger regions we cannot 
find extremal surfaces that stay outside the apparent horizon. This is not surprising 
since we expect based on earlier results that there will be surfaces that penetrate the 
apparent horizon of the black hole (cf., [24]). However, one of the disadvantages of our 
numerical scheme is that we are unable to explore this interesting regime due to the 
fact that the spacetime inside the apparent horizon has been excised. As explained 
in [4], this was to avoid complications with having caustics in the coordinate chart. 
Analysis of entanglement entropy however does require us to have the complete bulk 
geometry. 
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Fig. 16: Evolution of the extrennal surfaces for a strip of width a = 0.05 with driving parameters 
^{P = 10,^ = 20) = 200 (phase III; unbounded amplification). Conventions described in 
Fig. 14 apply. 

4.2 The evolution of entanglement 

We now turn to the evolution of the entanglement entropy; the results are presented 
in Fig 18 for the regulated quantity AS^ as introduced in (4.3). 

In the dissipation-dominated regime (phase 1), the entanglement entropy grad¬ 
ually increases, though in each cycle of forcing there is a time period for which the 
growth is negligible. We expect this feature is simply a consequence of the entangle¬ 
ment entropy tracking the thermal entropy. Even though we are not quite probing 
the full thermal contribution with the relatively small regions A, it bears to reason 
that the variation of the geometry is more or less equitable on all radial scales. This 
appears consistent with other probes of this phase. As we discussed in §3.1 the weak 
driving allows the system to efficiently dissipate the energy induced by the source 
and the conductivity a{t) was purely imaginary Basically the dominant effect here 
is the growth of the black hole horizon due to the driving and this in turn imprints 
itself into the growth of AS^ seen in Fig. 18a. 

On the other hand when we reach phase Ila (tilted regime) by way of small 
amplitudes, we start to see dehnite oscillatory evolution of AS^. In each oscillatory 
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Fig. 17: Evolution of the extrennal surfaces for a strip of width a = 0.01 with driving parameters 
^{P = 0.1,^ = 20) = 2 (phase Mb; dynamical crossover). Conventions are as described in 
Fig. 14. 


period we see a local reduction in AS^. On the other hand the temporal radial depth 
attained by the extremal surface as measured by u* is almost similar to that in phase 
I by juxtaposing the behaviour in Fig. 14 and Fig. 15. In phase Ila however, our 
extremal surfaces are closer to the boundary in contrast to phase I. We conjecture 
that the origin of the reduction in the AS^ is associated with the sharp oscillations 
in t* or equivalently i*. These imprint themselves into the actual value of the area 
despite the surface not getting too far into the bulk (which is possible since even the 
asymptotics of the geometry is sensitive to the driving, cf., (A.17)). The onset of 
non-monotone growth of AS^ in Fig. 18b characterizes the departure from the linear 
regime to the non-linear domain in line with the behavior of the phase portrait which 
in turns modifies the conductivity (which picks up a real part > 0 in phase Ila). 

The temporal change of AS^ is much more pronounced in the non-linear regime. 
In the unbounded amplification phase III (see Fig. 18c) and the dynamical crossover 
wobbly phase Ilb (see Fig. 18d), the AS^ appears to track the time-coordinate of 
the cap-off point i* quite efficiently. Indeed here we expect the non-linearities of the 
system to be the dominant effect. We know that the black hole grows quite rapidly in 
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(c) Phase III: ^{A = 20,P = 10) = 200 
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(d) Phase lib: ^{A = 20,P = 0.1) = 2 


Fig. 18: The evolution of the regularized entanglement entropy, AS^ defined in Eq. (4.3), for 
the four phases for a radial cutoff of Uj[ = 10“^. The strip widths are a — 0.05 for panels (a), 
(b), (c), and a — 0.01 for panel (d). 


response to the energy injected into the system at the boundary from our discussion in 
§3.2 and §3.3. The behaviour in phase III is smooth with large amplitude oscillations, 
which qualitatively track quite well the behaviour of i*. The dynamical crossover 
wobbly phase (phase Ilb) exhibits a lot more drastic behaviour. We encounter for 
the first time a jumps in the family of extremal surface that minimize the area 
(satisfying the boundary conditions and the homology constraint). These jumps 
translate into continuous but non-differentiable kinks in AS^ visible in Fig. 18d. We 
again note that the radial position of the cap-off point of the extremal surface behaves 
much more smoothly and the glitches appear in i*. Furthermore, the growth of the 
entanglement itself is rather steep as we see about an order of magnitude difference 
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Fig. 19: The evolution of the regularized entanglennent entropy, AS^ defined in Eq. (4.3), 
against the nornnalized entropy of the black hole, s/sq = s/s{t = 0), for the four phases for 
a radial cutoff of = 10“^. The strip widths are a = 0.05 for panels (a), (b), (c), and 
a = 0.01 for panel (d). We include the Spearnnan and Pearson rank coefficients, — I <Ps<l 
and — 1 < pp < 1 respectively, for each plot to demonstrate the linearity of the correlation 
between the entanglement entropy and the thermal entropy (see text for explanation). 


in AS^ between the low amplitude and high amplitude regimes. 

It is interesting to contrast the change of entanglement entropy with the change in 
the thermal entropy to see how the two are correlated. As we have argued above, the 
fact that we have an ever increasing thermal entropy (the bulk black hole is constantly 
growing) implies that even for small sub-systems we will quickly see overwhelming 
thermal contribution. We display in Fig. 19 the functional dependence of AS^ on 
the (normalized) instantaneous thermal entropy s{t)/s{t = 0). 
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It is immediately apparently by eyeballing the plots that there appears to be near¬ 
perfect correlation in three phases with Fig. 19c corresponding to phase III being the 
only outlier. To get a quantitative feeling for the correlation we have also indicated 
the Pearson correlation coefficient Pp as weli as the Spearman rank coefficient ps- 
These are statisticai markers for measuring correlations between two sets of data 
and are defined to take values in the interval [—1,1]. The values Ps,Pp = 0,±1 
signify zero, perfect positive and perfect negative correlation respectively. While the 
Spearman coefficient indicates that the observables in question are monotonically 
related, the Pearson coefficient provides an accurate measure of linear correlation. 
Indeed from the results quoted in Fig. 19 we see that AS^(s) is a linear function to 
a very good approximation in phases I, Ila and Ilb. It is curious that the linearity is 
respected even in the presence of the glitches in the growth of entanglement entropy 
(we do not see any drastic behaviour in the area of the apparent horizon). The 
unbounded amplification phase III clearly demonstrates the effects of non-linearities 
by decorrelating AS^ and s{t). 

5 Discussion 

The non-equilibrium dynamics of strongly coupled field theories is amenable to de¬ 
tailed quantitative exploration using the AdS/CFT correspondence. We have ex¬ 
ploited this set-up to study the behaviour when a homogeneous thermal plasma is 
driven away from equilibrium by a periodically sourcing a relevant (composite) scalar 
operator. The resulting dynamics exhibits a rather rich phase structure illustrated 
in Fig. 1. 

We identified four distinct phases, characterizing them in terms of the frequency 
and amplitude of the external driving force. Of these the dissipation dominated phase 
I is perhaps most intuitive for here the weakness of the driving, allows the system 
to to catch up with the driving. This is clearly visible in the various observables we 
studied; the complex conductivity of the response is purely real owing to the phase 
lag between the source and response and the evolution of entanglement is pretty 
quiescent. 

There is more structure when we ramp up either the period of driving, or the 
amplitude, for now the system departs quite rapidly away from equilibrium. The 
response therefore is more pronounced; we see more in phase response and greater 
temporal oscillations. In phases Ila to lib there emerges a non-vanishing imaginary 
part to the conductivity, which in fact appears to capture the entire response for 
high values of the period and amplitude. We also notice that there are significant 
fluctuations in the energy density and the entanglement entropy and furthermore, 
the entropy density grows rather rapidly in this regime. Perhaps most intriguing 
is the unbounded amplification of phase III, where we see sharp fluctuations and a 
highly non-linear response. We argue that this response appears to be not captured 
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by polynomial self-interactions of the composite operator; the intricate dynamics of 
gravity in AdS appears to induce effective non-polynomial couplings in the effective 
action for the operator O we use to perturb the system away from equilibrium. We 
believe this fact is significant and should be taken into account when attempting 
to construct effective models distilling the effects of gravitational interactions for 
strongly coupled systems . 

While our focus has been on computing the simplest set of observables, essentially 
one-point functions and entanglement entropy for small sub-systems, the power of 
holography is that we can do much more. In time independent equilibrium scenarios 
it is straightforward to use the holographic map to compute correlation functions (at 
least two point functions). In the genuine non-equilibrium scenarios as those we have 
focused on the technology for computing such observables, whilst present [45, 45] is 
still a bit cumbersome to work with (at least numerically). It would be interesting 
to develop these techniques further perhaps taking inspiration from the analytical 
models of [54, 55]. This would allow us with a direct probe of fluctuations in the 
plasma, which can be contrasted with the dissipation in the system, the latter being 
measured by the entropy production through the growth of the horizon. 

Likewise our exploration of the behaviour of entanglement entropy has been re¬ 
stricted to analysis of small sub-systems for pragmatic reasons. While the sub-system 
under consideration was chosen to have fixed size, the fact that we are continuously 
driving the system leads to an ever increasing thermal contribution to the entangle¬ 
ment. Geometrically this is easy to understand since the horizon for our bulk solution 
is ever growing (as we have indicated that both the event and apparent horizons are 
required to be monotonic in our set-up) and reaches out towards the boundary in 
the course of the evolution. As a result, the local thermal scale can overwhelm the 
relative smallness of the sub-region we choose. To have precise mapping of the en¬ 
tanglement structure we need to be able to ascertain the true minimum of the area 
functional in such scenarios bearing in mind that the extremal surface can (and often 
does) penetrate various horizons. A significant obstacle in ascertaining this is the 
fact that the characteristic method for solving Einstein’s equations developed in [4] 
excises the region of the spacetime behind the apparent horizon. While this is a 
technical obstacle, overcoming it would not only enable us to probe the interior of a 
highly non-equilibrium black hole using holographic entanglement, but it could also 
allow us to explore other interesting scenarios such as the effect of perturbing the 
ground state of the system by external sources. 
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A Holographic Renormalization 

We collect here some salient results for the computation of physical field theory 
quantities using standard holographic techniques. 


A.l Scalar deformations 

The bulk action (2.5) should be supplemented by boundary counter-terms to ensure 
that (a) the bulk equations of motion follow from a consistent variational principle 
and (b) the on-shell action evaluated on the solutions is finite. 

In standard Poincare-AdSd+i 


ds^ 


dx^ dx'' -F — 


rini, dx^ dx^ + dz^ 


(A.l) 


the scalar field behaves asymptotically as 


(f){z,x) ^ (pod- 4)1 (A.2) 


We will work with standard quantization (Dirichlet boundary conditions) for the 
scalar field, which involves treating the mode that fall-off as as the source for 
the scalar held. 

In the presence of the source we let the metric to take the FG form, 

= 9nv{x,z)dxi^ dx'^ 

where gi^i,{z,x) = -|- 0{z). If necessary we will denote by % the induced metric 

on the surface z = z^ which differs from the boundary metric by a conformal trans¬ 
formation by z^. We will ignore this issue for most part and write the counter-terms 
in terms of below for simplicity. 
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With these conventions we find the following boundary counter-terms: 


S, 


bdy 


= ( 2W - 2 (d - 1 ) - 

IfivrGjv 7 V / y V ’ d-2 


m 


1 2 1 
--A_0 + 2(2A-d-2) 

We are using conventional AdS/CFT definitions: 


[(90)2 +Cl 


(A.4) 


d 


d2 


d 


A+ = -±W—-hm2 = -±i/ 


(A.5) 


Our interest concerns conformally coupled scalar field which has a mass in AdS 
units given by 


ml = 


d? - 1 


A+ = 


(i zb 1 


4 ^2 

To compute the boundary energy momentum tensor we vary 


rpilU _ 


2 -\- Sjjfly) 


(A.6) 


(A.7) 


where we should take care to include the appropriate radial dependence in the defi¬ 
nition of 7 ^ 1 ,. 

Lets split the contribution from the graviton and the scalar and write 




(A.8) 


where the split is determined by the requirement that oc 0. Then the two pieces 
can be computed efficiently as follows: 


IfiTTG^ C=7 H 


flU 


d'^+^x ^f^{R + d{d-l)) 


f 


+ / d^x 2 iF - 2 (d - 1) - 


1 


d-2 




which one can show evaluates to a nice covariant expression: 


_ _ 

^ “ IfiTTG 


N 


+ {d- 1) 7^^ - - -^^7^'^ 

d 2 \ 2 


(A^9) 


(A.IO) 


where is the location of the cut-off surface. 

The scalar contribution can be evaluated by using the fact that we are interested 
in the boundary variations to obtain: 


rp/diy _ 


1 


IfivrGjv ^7 


IIU L' 
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2zf 


(A.ll) 
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where • • • indicate the contribution from the higher order counter-terms and we have 
put back the powers of now. The details now depend on the asymptotic expansion 
of (p. For general A we have to worry about the fact that the Taylor series solution 
in the neighbourhood of z ~ 0 looks like 

0(z, x) = 00 + Oi(0o) z^-~^ -\ -1-01 Z^+ -\ - (A.12) 


and we need to know the various intermediate pieces to complete the analysis. The 
case we are interested in is rather special, where there are no powers of z in the 
Taylor expansion between the source and the vev, so let us simply record the result 
for this case for now leaving a more general analysis for later. 

Before proceeding though, let us note that we can express (A.11) covariantly as 
follows (r = z~^)\ 


rjifiy _ 


IbvrGjv hi 


flU 




(A,13) 


where h is the unit normal perpendicular to the cut-off surface. 


A.2 Specializing to A_ — A+ < 2 (A+ > d — 2) 

In this case the asymptotic expansion belongs to the special kind where 

(l>{z, x) = 00 z^~ 01 z^"*" -!-••• (A.14) 

where we are allowed to use the fact that 0o z^“ is the beginning of an independent 
Taylor series where the powers of z change by 2 units (use the fact that the Lagrangian 
has 0 —7- —0 symmetry). This corresponds to the case we are interested where 
A+ = 2, A_ = 1 in d = 3. 

In this circumstance we can simply use the terms written explicitly in (A. 11) to 
obtain 

Then we find 

-Kj^^ + {d- 1) 7^^ + ^ (2A+ - d) 0o 0i 7 m^^^ (A-16) 

A.3 rh = —2 in d = 3 


Now, we can get the final answer for the case of interest either by working with the 
Fefferman-Graham expansion in which case we need to know that 


Qnuiz, x)dx^ dx" 


+ “ 4 ^0 


asz 


dt^ 


~ 2 (^3 + 00 0l) + 


{dx^ -|- dy^) (A.17) 
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The metric fall-offs allow us to compute the pieces in the boundary stress tensor 
directly since the term above is the correct answer. 

Using this or directly computing from the CY-ansatz (A. 13) we claim to obtain 
(rescaled the result by a factor of 3/2). 

= diag|2o300 01,- 03 ,-osj (A.18) 

We can check that this satisfies the Ward identities: 

TH = (j)0(P, = JO2, = -2 03 - 01 00 - 00 01 = -01 = (A. 19) 

where we used the boundary conservation law derived from the solution 03 = — ^ 0 o 0 i • 

B Extremal surfaces and entanglement for strips 

In this appendix we describe our methodology for finding extremal surfaces relevant 
for the computation of entanglement entropy. For simplicity we will focus on regions 
which exploit the symmetry of our set-up and consider A to be a strip extended along 
one of the translationally invariant directions, say y without any loss of generality, 
as in Eq. (4.1). We need a bulk codimension-2 surface that ends on the boundary of 
this region i.e., at x = ±0 (at the chosen instant of boundary time tj\). We describe 
our strategy for finding this surface and computing its (regulated) area below. 

B.l Determining extremal surfaces 

To find the extremal surface, we start by gauge fixing the reparameterization invari¬ 
ance on the surface. We take y to be one of the coordinates. Dimensionally reducing 
in this direction, we construct an effective action for a curve in the remaining di¬ 
rections and pick a proper-length parameter A as the second coordinate. Thus, the 
extremal surface Sji, is embedded in the bulk as 

Y^ = (t(A),r(A),x(A), 2 /). (B.l) 

We choose the proper-length parameter to ensure that \/det00/ = 1, which implies 
that the unregulated area of the extremal surface is given as 

Area(^^) Ly I dX\J dety^^^ , (^’^) 

Jsa 

in terms of parameter distance A^^ spanned by the curve and the IR regulator Ly. 
In practical terms we work with the effective Lagrangian 

C = [ 2 1' (r' — t' f) + (B.3) 
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where the metric functions p, f, x obtained by interpolation of our numerical 
solutions. This is a geodesic problem, with some non-minimal coupling from the 
dimensional reduction along the translationally invariant direction of the strip. In¬ 
stead if using the geodesic equations, we found it convenient to pass to a set of six 
first-order Hamilton-like equations by introducing Pt = f, = x\ and P+ = r' — pt' 
which are related to the conjugate momenta. The equations we solve are the above 
three and 

4P 

K = iW + P+)drP + Ptdtp) = 0 

2 p2 

p; = 2P^pe-^^drP - p2 (drf + 2fdrX + 2dtx) - -y ifdrP + dtp) = 0 (B.4) 

2 p2 

p; = 2P^pe-^^ {fdrP + dtp) + PtP+drf - {drP + pdrX) = 0 

We start from x = 0 in the bulk at some smooth cap-off point (x = 0,t*,r*) where 
t' = r' = 0.^^ and propagate out to the boundary. We evolve until a with a hxed 
UV cut-off at and regulate the final answer for the entanglement entropy by 
background subtraction (see below). 

In the main text we illustrate the temporal dependence of the extremal surfaces 
and for each of the four phases (I-IV) of Fig. 1 for fixed strip width a. Since 
we numerically control the data of the cap-off point we work iteratively: we start 
by fixing a suitable strip width a by tuning r* and t*, then we evolve the extremal 
surfaces by increasing t* and re-adjusting r* such that the strip width remains as a. 
We note that we assume that there are no discontinuities or multi-valuedness in the 
map from —)■ (a,r_ 4 ), which we believe makes sense for small strip widths.^® 

Finally, to work in a compact domain we choose u = 1/r e [0,1] which we will use 
to explain the properties of the extremal surfaces. 

B.2 Regulated entanglement entropies 

Since the extremal surfaces reach out all the way to the boundary, the proper area is 
divergent with the coefficient of the leading divergent term fixed by the area of the 
entangling surface dA. For a state of the CFT with vanishing sources for operators 
it is well known [52] that the entanglement entropy behaves as 

= + + (B.5) 

where is finite in the limit ^ 0. In normalizable states of the field theory 
is the universal contribution to entanglement which should be independent of the 

This cap-off point is not necessarily the deepest point in the bulk; for the examples shown in 
this paper it however does turn out to coincide. 

Such behaviour was noticed in extremal surface computation in global Vaidya-AdS by [53] . 
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cutoff value One natural way for us to extract this quantity is to measure the 

entanglement relative to the t = 0 thermal Schwarzschild state ASj,{t) = S_A{t) — 
Sjiit = 0), which can be extracted simply by vacuum subtraction. 

Usually, when we turn on sources for relevant operators, these can contribute 
additional divergences to the entanglement entropy [60]. In general in the presence 
of additional relevant scales one naively expects there to be logarithmically divergent 
terms polluting (B.5) and rendering vacuum subtraction meaningless. Fortuitously, 
this does no happen for the problem at hand. This can be extracted by examining 
the detailed discussion of [60], which we paraphrase below. 

There is however a quick argument for the absence of logarithmic terms which we 
now describe. For scalar operators in CFT^ with operator dimension I < A < I+ 1, 
as we have considered, it is well known in AdS/CFT that the corresponding bulk 
field has mass in the window where both asymptotic fall-offs are normalizable, i.e., 
m? G + 1) with the Breitenlohner-Freedman bound mass m%p = — ^ as 

usual.In this window note that A_—A_|_ < 2 and we have the Legendre transformed 
theory with an operator of dimension A_ by switching to alternate quantization [61]. 

Turning on a source for the faster-fall off mode A+ is equivalent, insofar as the 
leading back-reaction on the metric, to considering instead a state in the Legen¬ 
dre transformed theory where the alternate quantized operator with dimension A_ 
acquires a vacuum expectation value. However, since the divergence structure of 
entanglement is the same in all states of the field theory, and the conformal vacua of 
the two theories (standard and alternate quantization) coincide, it follows that the 
divergence structure of S_a should be unchanged from (B.5), even with J7(x) ^ 0. 
Our story is of course a special case with A+ = 2, A_ = 1 in d = 3. This observation 
is consistent with the results of [60] and the counter-terms used in [20]. 

To explicitly analyze the structure of the divergences in the entanglement en¬ 
tropy, let us consider the metric given in (A. 17). Since the details of the divergences 
are blind to the boundary spatio-temporal behaviour of the sources we will examine 
the somewhat simplified setting where 00 = const to glean the relevant information. 

With the time-translational symmetry restored by this choice, the Lagrangian 
for the extremal surface (which now is minimal) is simpler: 

^ \/gii{.z) + z'{xy (B.6) 

For the vacuum state of a CFT 3 with A being a circular disc would give the F-function 
[56, 57] (the latter defined as the logarithm of the partition function of the theory a three-sphere). 
In fact, this can be used to define a UV finite quantity without recourse to background subtraction: 
following [58, 59] we can just as well consider (i? ^ — l) 5'^, with R being the disc radius, as the 
measure of entanglement growth. 

Implicit in this statement is the fact that we are quantizing the scalar field with standard 
(Dirichlet) boundary conditions, so that the dimension of the dual operator is A = A+. 






where gii{z) is the spatial component of the metric in (A.17). This system has a 
conserved Hamiltonian, which we can exploit to write down an expression for the 
area directly. Introducing, z* which captures depth to which the minimal surface 
penetrates into the bulk, we have for the on-shell value of the area 


Area oc 






1 

2 


(B.7) 


Using the explicit form of gu^ the second term is at least near the boundary so we 
can forget about it. The first term is all that matters, so lets look at 


^ — 2 (<^3 + ^ + • • • (B.8) 

which has the divergence expected upon integration, but no further contribution 
of relevance in z ^ 0 limit. From the (/)q term we get a contribution to the finite part 
of the entanglement, and this is indeed the physically relevant answer. It should be 
clear from this discussion is not specific to the choice = — 2 in d = 3, but should 
hold for I < A+ < | + 1 as we argued abstractly above. 
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